clc;
clear all;

%%
load('section4_all.mat');
format short g
%% 
% unpack data needed
donation=d.donation; default=d.default;
relfreq=d.relfreq;

n=length(relfreq);
% calculate Delta: gain of deviating from default
Delta=Delta_fun(donation, default);

%%
% exclude AD group for analysis
Delta= Delta(default~=0);
donation= donation(default~=0);
default= default(default~=0);

%% Set seed for replication. 
% the global search algorithm uses a pseudorandom number generator to pick 
% multiple starting values to solve for the  minimum. This reduces the
% chance of ending up in a local minimum. 
% Because of the randomness, without setting the seed the algorithm returns parameters 
% that are nearly identical but not exactly the same as the results in the paper. 
% By setting the seed, the exact same results can be obtained. 

rng(432123486);
%% Use global search with fmincon
% constraint optimization with:
%- lambda1 between 0,1
%- lambda2 and lambda3 nonnegative
% global search starts a local server from multiple starting points
% no rescaling, set c1 and c2=1
likelihood=@(lambda)totallikelihood(lambda, donation, Delta, default, relfreq, 1, 1);
opts=optimoptions(@fmincon,'Algorithm','interior-point', 'FunctionTolerance',1e-12, 'StepTolerance',1e-8,'AlwaysHonorConstraints','bounds');
problem = createOptimProblem('fmincon','x0',[0.8; 0.0005; 0.005], 'objective', likelihood, ...
  'lb', [0 ; 0; 0], 'ub',[1; Inf; Inf],'options', opts);
gs = GlobalSearch;
[x,FVAL,exitflag] = run(gs,problem)


display('estimated values for lambda');
disp(x');
%% save estimated lambdas in mat file
save('lambdas_hat.mat', 'x');


%% When no access to the global search toolbox, fmincon can be used
% however note that because different starting values are used, the results
% may differ. 
%{
% fmincon
% Goal: recover parameters used for simulation 
% lower bound
lb=[0 , 0, 0]; % constraint optimization, lambda1>=0, lambda2>=0
% upper bound
ub= [1, Inf, Inf]; % lambda1 constraint: between 0 and 1
lambda0=[0.8, 0.0005, 0.0005]; % starting values 

% no rescaling 
c1=1; c2=1;

options = optimset('Display','iter','GradConstr','off','GradObj','off','MaxIter',10000,'DerivativeCheck','on','MaxFunEvals',2000);
options1 = optimset(options,'Algorithm','interior-point', 'TolFun',1e-12, 'TolX',1e-8,'AlwaysHonorConstraints','bounds');

[x1,FVAL1,exitflag]= ...
    fmincon(@(lambda)totallikelihood(lambda, donation, Delta, default, relfreq, c1, c2), lambda0, ...
    [], [], [], [], lb, ub, [], options1);

%}

%% Numerically approximate BHHH Hessian: 
% Use variant of BHHH algorith to numerically approximate hessian, 
% since approximated Hessian from fmincon is not accurate. 
% BHHH hessian is positive definite by construction. 
% Avoid problem of invertibility of hessian when parameter estimates are
% small.
% Bt= sum (Sn(lambdat)Sn(lambdat)')/N average outerproduct of the sample
% Sn is the score, Sn(lambdat)=dLoglik(lambdat)/dlambdat
% the average outproduct of the sample is related to the covariance matrix
% of scores in the sample. The score is an observation sprecific gradient and
% Bt is the average outerproduct of these observation specific gradients
% at the parameters that maximized the likelihood function, the average
% score is zero. then the outerproduct of the scores, Bt, becomes the
% variance of the scores. 
lam_true= x
% numerically approximate score
score1= numscore(lam_true,  donation, Delta, default, relfreq, 1, 1,  1, 0.000001);
score2=numscore(lam_true,  donation, Delta, default, relfreq, 1, 1,  2, 0.000001);
score3=numscore(lam_true,  donation, Delta, default, relfreq,1, 1,  3, 0.000001);
Sn= [score1' ; score2' ; score3'];

% BHHH Hessian is outerproduct of scores
BHHH= (Sn*Sn');
SD1_bhhh = 1/(sqrt(BHHH(1,1)));
SD2_bhhh= 1/(sqrt(BHHH(2,2)));
SD3_bhhh=1/(sqrt(BHHH(3,3)));
SD_bhhh= [SD1_bhhh, SD2_bhhh, SD3_bhhh];


%% pack in data structure 

lambda= x; sd_bhhh=SD_bhhh'; likelihood=[FVAL, NaN, NaN]';

%%
T=table(lambda, sd_bhhh, likelihood);

writetable(T,'results_all.xls') 

%% for results with rescaled parameters, unblock the following code
%{ 
%% rescale lambda2 and lambda3, set c1 and c2=100
clc;
clear all;
load('section4_all.mat');
format short g
%%
% unpack data needed
donation=d.donation; default=d.default;
relfreq=d.relfreq;

n=length(relfreq);
% calculate Delta: gain from deviating
Delta=Delta_fun(donation, default);

%%
% exclude AD group for analysis
Delta= Delta(default~=0);
donation= donation(default~=0);
default= default(default~=0);

%% global search

likelihood=@(lambda)totallikelihood(lambda, donation, Delta, default, relfreq, 100, 100);
opts=optimoptions(@fmincon,'Algorithm','interior-point', 'FunctionTolerance',1e-12, 'StepTolerance',1e-8,'AlwaysHonorConstraints','bounds');
problem = createOptimProblem('fmincon','x0',[0.8; 0.5; 0.5], 'objective', likelihood, ...
  'lb', [0 ; 0; 0], 'ub',[1; Inf; Inf],'options', opts);
gs = GlobalSearch;
[x,FVAL,exitflag] = run(gs,problem)


display('estimated values for lambda');
disp(x');

%%
%{
% fmincon
% Goal: recover parameters used for simulation 
% lower bound
lb=[0 , 0, 0]; % constraint optimization, lambda1>=0, lambda2>=0
% upper bound
ub= [1, Inf, Inf]; % lambda1 constraint: between 0 and 1
lambda0=[0.8, 0.5, 0.5]; % starting values 

% rescaling 
c1=100; c2=100;

options = optimset('Display','iter','GradConstr','off','GradObj','off','MaxIter',10000,'DerivativeCheck','on','MaxFunEvals',2000);
options1 = optimset(options,'Algorithm','interior-point', 'TolFun',1e-12, 'TolX',1e-8,'AlwaysHonorConstraints','bounds', 'MaxFunEvals',2000);

[x1,FVAL1,exitflag]= ...
    fmincon(@(lambda)totallikelihood(lambda, donation, Delta, default, relfreq, c1, c2), lambda0, ...
    [], [], [], [], lb, ub, [], options1);

%}

%% Numerically approximating BHHH Hessian: 

lam_true= x
score1= numscore(lam_true,  donation, Delta, default, relfreq, 100, 100, 1, 0.000001);
score2=numscore(lam_true,  donation, Delta, default, relfreq, 100, 100,  2, 0.000001);
score3=numscore(lam_true,  donation, Delta, default, relfreq, 100, 100,  3, 0.000001);
Sn= [score1' ; score2' ; score3'];

BHHH= (Sn*Sn');
SD1_bhhh = 1/(sqrt(BHHH(1,1)));
SD2_bhhh= 1/(sqrt(BHHH(2,2)));
SD3_bhhh=1/(sqrt(BHHH(3,3)));
SD_bhhh= [SD1_bhhh, SD2_bhhh, SD3_bhhh];


%% pack in data structure 

lambda= x; sd_bhhh=SD_bhhh'; likelihood=[FVAL, NaN, NaN]';

%%
T=table(lambda, sd_bhhh, likelihood);

writetable(T,'results_all_rescaled.xls') 
%}
